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1. Introduction 

The problem of determining the orbit of a space object from measurements based on one 
pass through the field of view of a radar is not a new one. Extensive research in this area has 
been carried out in the USA and Russia since the late 50s when these countries started the 
development of ballistic missile defense (BMD) and Early Warning systems. In Russia these 
investigations got additional stimulation in the early 60s after the decision to create a Space 
Surveillance System, whose primary task would be the maintenance of the satellite catalog. 
These problems were the focus of research interest until the middle 70s when the appropriate 
techniques and software were implemented for all radars. Then for more than 20 years no new 
research papers appeared on this subject. This produced an impression that all the problems 
of track determination based on one pass had been solved and there was no need for further 
research. 

In the late 90s interest in this problem arose again in relation to the following. It was 
estimated that there would be greater than 100,000 objects with size greater than 1-2 cm 
and collision of an operational spacecraft with any of these objects could have catastrophic 
results. Thus, for prevention of hazardous approaches and collisions with valuable spacecraft 
the existing satellite catalog should be extended by at least an order of magnitude This is 
a very difficult scientific and engineering task. One of the issues is the development of data 
fusion procedures and the software capable of maintaining such a huge catalog in near real 
time. The number of daily processed measurements (of all types, radar and optical) for such 
a system may constitute millions, thus increasing the number of measurements by at least 
an order of magnitude. Since we will have ten times more satellites and measurements the 
computer effort required for the correlation of measurements will be two orders of magnitude 
greater. This could create significant problems for processing data close to real time even for 
modern computers. Preliminary ’’compression” of data for one pass through the field of view 
of a sensor can significantly reduce the requirements to computers and data communication. 
This compression will occur when all the single measurements of the sensor are replaced by 
the orbit determined on their basis. The single measurement here means the radar parameters 
(range, azimuth, elevation, and in some cases range rate) measured by a single pulse. 

Two types of techniques have been traditionally used for processing single measurements; 
recursive and joint or batch processing. Recursive procedures convenient for real time pro- 
cessing usually are based on Kalman’s recursive filter [1], Less convenient joint processing 
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techniques basically use the least squares [2] or least modules methods [3]. When the single 
measurement errors are time correlated, and when the statistical characteristics of the sin- 
gle measurements are not known completely these techniques do not provide a guaranteed 
evaluation of the errors of the generated estimates. 

This limitation can be avoided when we use the joint method with the guarantee approach 
when the guarantee ranges of orbital parameters are obtained on the basis of guarantee ranges 
of the parameters of single measurements [4], The guarantee approach has one more remark- 
able feature. With certain limitations on the distribution of the errors of the measurements 
and with a large enough number of these measurements this approach leads to a much more 
accurate estimate than the traditional techniques mentioned above. 

The general procedure based on the guarantee approach is a linear programming procedure 
with the amount of computation significantly greater than the least squares procedure. This 
resulted in a lack of interest to such procedures in 60-70s. However, now the situation is 
different, the capacity of the computers is significantly greater. Modern sensors have small 
and rather stable errors. Thus we have the reasons to look again at this promising method. 

The current paper presents the comparative analysis of the accuracy characteristics of 
different procedures using mathematical simulation with the following background data: 


• Limitations for the composition of the measured coordinates: 

local spherical coordinates : d, a, (3 (1) 

• Limitations of the field of view of the radar: 

150 km < d < 7000 km 0 < f3 < 50° (2) 

• Limitations for the accuracy of the not abnormal single measurements: 

0.01 km < a d < 0.05 km 0.01° < cr a , o g < 0.03° (3) 

• Limitations for the time interval between neighboring measurements: 

1 s < a < 10 s (4) 

• Limitation for the tracking time: 

At < 300 s (5) 

• Limitations for the observed satellites: 

no active operations and small atmospheric drag (6) 

• Limitations on satellite parameters: 

i > 30° h p > 100 km e < 0.8 (7) 
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2. The Joint (non-recursive) Algorithm Based on the Least Squares Method 


The least squares method search for the mind/ (a) = d/ (a rmn ) of the function dRa) of the 

a 

form 

n / i i i \ 

dRa) = VW — (4 - 4(a)) 2 + — (a fc - a fc (a)) 2 + — (4 - /4(a)) 2 , (8) 

W < J 

where 


• x fc = (4, cx k , Pk) ~ k-th measurement (k = 1,2, ...,n); 

• t k - time of the k-th measurement; 

• &d k , cr ak , a 9k - RMS of the errors 4, a fc , 4; 

• h fc (a) = (4(a), Ofc(a), /4( a )) - values of the parameters of the k-th measurement, calculated 
using the vector a = a (t) of the orbital parameters referred to certain time t. 


For the case of Gaussian non-correlated errors of the measurements the estimate a mm pro- 
vides a maximum for the probability density function p (x 1; x 2 , ..., x n ) , thus it is the maximum 
likelihood estimate having the feature of asympthotical efficiency [5]. The covariation matrix 
of the errors K of the estimate a mm in this case can be calculated using the formula 


K 


0.5- 



(9) 


The estimate a min retains the asymptotic efficiency feature for the case of low level 
correlation between measurements (correlation interval significantly smaller than the interval 
of observations). If all the errors of the measurements are Gaussian and all the functional 
relationships x fc (a) are linear, then the estimate dL m .i n is unbiased and has the least RMS errors 
within the class of unbiased estimates having the shape of linear functions of measurements. 

These remarkable features of the least squares technique make it attractive for applications 
in practice. However, implementation of this method sometimes faces several problems. The 
major ones are: 

- choosing the vector a and the method for its propagation for time interval r; 

- calculation of the initial approximaiton a 0 , for minimization of T(a); 

- selection of technique for reaching the minimum of dRa). 

For the vector of estimated parameters a it is convenient to use the state vector 
(x, y, z, x, y, z ) for the time period in the local rectangular coordinate frame (LRCF) related 
to geographical directions of the sensor. In this coordinate frame the basic plane is the plane 
of the local horizon, the x axis is in the plane of the local horizon and is directed to the west, 
the y axis is directed normal to the plane of local horizon upwards and the z axis is in the 
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plane of local horizon and is directed to the north. Coordinates x, y, z can be obtained 
from local spherical coordinates d, a, (3 using the formulas 

x — —d sin a cos [3 y — dsm/3 z = dcosacos/3 (10) 


The equations of motion in LRCF with accuracy sufficient for solving the task mentioned 
in the introduction have the form [2]: 


d 



— 2cuxd — ci;x(a;xr) 


3 J 2 yR 2 e 

2 r 5 


fr + 2(r,cJo)cUo 



( 11 ) 


where 


• d = (x, y, z)', d = (x, y, z)' - position and velocity vectors in local coordinate frame; 

• r = {x—b x , y—by, z—b z )' - position vector of an object with respect to the center of the 
Earth in local rectangular coordinate frame; 

• r — \J (x— b x ) 2 +(y— b y ) 2 +(z— b z ) 2 - length of the vector r; 

• b = (b x , by, b~y - coordinates of Earth’s center in local rectangular coordinate frame - con- 
stants, which values are determined in the following way using the coordinates X g0 , Y g0 , Z g0 
of the radar locaiton (the origin of local rectangular frame) in the Greenwich rectangular 
frame: 


Jxio+Yl* r 

= \l X lo+ Y io+ Z & z = Z g0 /r f l = rjr 

(12) 

R e -(l—a-z 2 ) 

b x = 0 b y = —r b z = 2-a-R-z-fi 

(13) 


a = 0.0033528107 - is the constant characterizing Earth oblateness 


• u — ctvcuo, where u) e - angular velocity of Earth’s rotation, dy, - unit vector of the 
angular velocity of Earth’s rotation in local rectangular coordinate frame, which components 
ojy x , oj^y, u)y z are the constants calculated using the formulas: 

^ = 0 u oy = z + n-bz/r uj oz = iq — z-b z /r (14) 


• y = 398600.44 km 3 /c 2 - gravitational constant; 

• J 2 = 0.001083 - coefficient of the second harmonic of the Earth’s potential expansion; 

• R e = 6378.137 km - equatorial radius of the Earth; 

• (a, b), axb - scalar and vector products of vectors a and b; 

• the vectors are denoted in small letters, matrices - in bold ones; 

• vector with Latin notation are in bold letters, with Greek one - have superscript; 
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• sign ' after vector or matrix means transposition; 

For solving these equations (propagation of the state vector to the given time r) we use 
the numerical 4th order Runge-Kutta method. 

For solving the system of differential equations a (t) = F(a(f)) c a(t=0) = a 0 this 
method uses a stepwise process a i+1 = a, + Aa*, where a i+1 - is the result of propagation of 
the parameters a 0 for the time T 1 +r 2 +...+Tj until we reach the time t — r. For the i-th 
step the increment Aa* is calculated by the formula: 

Aa,; = (kf> + 2k f + 2k f + kf)/6, (15) 

where 

kp=Ti-F(ai) k^=r r F( ai +0.5k^) k^=r r F( ai +0.5k^) kf=7yF( ai +k<*>) (16) 


Selection of the technique for the minimization of T(a) that will have a guaranteed and 
quick convergence is not an easy task, ffowever, the condition (2) generates the situation 
when during the tracking interval the major effect of updating by least squares is the improved 
accuracy of the velocity components u=(i, ij, z ) of the state vector a. The position 
parameters u=(x, y, z ) of the vector a in fact are not updated. Thus we suggest for 
the minimization of T(a) = T(u, u) = -0(u) we use only the vector u. Under the limits 
(3) of the measurement errors and the technique for calculation of the initial approximation 
selected below the relationship x(u) is close to linear. Thus the iteration process converges 
very quickly. Eventually for all the cases the search for the minimum u m ; n of the function 
■^(u) requires only one iteration. 

For the initial approximation a 0 We can use the estimate a n , obtained using measure- 
ments Xj, x 2 , ...,x n by the recursive algorithm described in the next section. 


For the point of the minimum u mtn of the function ^(u) for the case when it is obtained 
in one iteration we suggest using Newton’s formula 


U min 'hi 




( 17 ) 


where the Erst and the second derivatives of the function ^(u) with respect to parameters 
x, y , z are calculated using finite differences. 


3. Recursive Algorithm 

We know that the estimate a min of the least squares method has a recursive structure: 
the estimate of the parameters based on k measurements can be written as a function of 
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the same estimate for k — 1 measurements and the k-th measurement. This is true when all 
the functions h fe (a) are linear. For non-linear problems including our case this is sufficiently 
accurate when the errors of the measurements are relatively small and a satisfactory initial 
approximation for orbital parameters is available. 

The recursive feature of the least squares estimate can be extended to the more general 
case [1], when the parameters a fc =a (t k ) Of the system for the time t k (in our case the orbital 
parameters for the time t k ) in addition to deterministic component have a random one which 
has Markov’s feature (the ’’future” for the fixed ’’past” depends only on the ’’present”), i.e. 

P (a fc |a fc _i, Ufc- 2 , •••) = P ( 3 * | fJ-fc—i ) (18) 

where p(a k \....) - conditional probability density a fe . 

The condition (18) is satisfied if 

a k = i k (a fc _i) + (19) 

where y fc - ’’noise of the system” - time independent random perturbations with zero mean 
and covariation matrices T fc . 

Regarding the problem under consideration these perturbations may be related to the errors 
of the propagation model (in particular the atmospheric density model) and the calculations 
(e.g. roundoff) errors. 

For (18) and independent measurement errors of the aposteriori (after acquisition of the 
measurement x fc ) probability density of the parameters a fc we have the following recursive 
formula 


p(a»|x„ — p (a fc |x fcl a,_,) =p(x,|a 4 ).j>(a,|, a,_ 1 )/p(x,), 


( 20 ) 


From this relationship under the assumption (19) we can derive the approximate recurrent 
formulas for the first and second moments a ki P fc of the aposteriori distribution a fe 


where 


P,: 1 = P*iLi + H 1 R * H a p fc % = Pfcifc-i^fc|fc-i + 


a ( A s p _df k (&k-i) tt _^hfc(^fc|fc-i) p _p p p/ | p 

— L k\ <± k-l) r k — 7^ 7^ A k\k-l — r fcA k-l r k 




( 21 ) 

( 22 ) 


Equations (21), (22) are exact if the functions h fe (a fc ), f fc (a fc _,), relating the measured pa- 
rameters to the orbital parameters for the times t k and t k _ 1 , are linear, and the distributions 
of the measurement and dynamic model errors are normal. 

Equations (21) and (22) present a non-linear generalization of the known recursive 
Kalman’s filter. For our case when the measured parameter is a scalar u* 7 the matrix 

*We can consequently update the previous estimate using the current measurements of range 
(u = d ), azimuth (u = a) and elevation (u = /3), since the measurements d, a , (3 are not correlated 
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H fr becomes a vector-column h fc and (21) are transferred to the form 

Wk = w k /( 1 + w k 

Pfc = P k\k~i — ^fc(P fc|fc-ih fc )(P fc|fc_].h fc ) (23) 

3-k — aib|fc-i + WkPfch k (u k — h fc (a fc | fc _ 1 )) 

where w k = 1 /(J 2 Uk - is the weight characteristic of measurement u k . 

The choice of the specific algorithm for processing the measurements assumes that the 
following should be determined: 

• parameters u k —u(t k ) of the measurements for the time t kr 

• parameters a fc =a (t k ) of the orbit for the time t ki 

• the operator f k (a fe _ 1 ) for propagating the orbital parameters from the time t k _ 1 to the 
time t k , 

• operator F fe (a fc _j) for the propagation of the variations of orbital parameters from time 
t k _ i to the time t k , 

• the technique for generating the initial approximaiton a 0 , P 0 and t 0 for calculating using 
the recurrent relationships (16), 

• the technique for calculating the noise matrix of the system I\. 

The coordinate frames of the parameters u k , a fe and the propagation operator f k (a fc _ 1 ) 
are the same as selected for the least squares method. 

For the operator F fc (a fc _i) we will take the matrix F fe with dimensions 6x6 with the 
following non-zero elements: fi,i= 1 for i = 1, 2, 6, /i,i +3 =r fc fori = 1,2,3, where 

Tfc = tk~t k _ i . 

1 0 0 T n 0 0 

0 1 0 0 T. n 0 

0 0 1 0 0 T n 

0 0 0 1 0 0 

0 0 0 0 1 0 

0 0 0 0 0 1 

The initial approximation a 0 = (x 0 , y 0 , z 0 , x 0 , y 0 , z 0 )' = (d 0 , d 0 ) for the vector of orbital 
parameters for the time t 0 and the matrix P 0 of its errors are determined from the first two 
measurements x, = (d 1; a u /3J and x 2 = (d 2 , a 2 , /3 2 ), acquired for the times t 1 and t 2 . 
Regarding the limitations (3) and (4) the following procedure will be correct: 

• The measurements x, , x 2 are transferred to d, , d 2 using formulas (10). 

• Then we calculate the state vector (d 0 , d 0 ), consistent with the limitations (3) for the 
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measurement errors of the angular coordinates: 


d o =0.5(d 1 +d 2 )-r 2 d o /8 d 0 =(d 2 -d 1 )/r t 0 = 0.5(t 1 +t 2 ) 


(24) 


where r = (t 2 — H), a d 0 is determined using (11) from the state vector without compen- 
sating term r 0 2 d 0 /8. 

In the calculated state vector (d 0 , d 0 ), the parameters d 0 , d 0 are changed to more 

accurate values d 0 =\fi o, d 0 =£ 0 /2d 0 , adequate to the limitations (3) for the errors of range 
measurements: 

lo=0.5(e i +6)-r 2 e o /8 l=(6-e i )/r-r 2 e ) /24 (25) 

where £=d 2 , £ 0 =2(d 0 , d 0 )+2(d 0 , do), £o 3)= 6(do, d 0 )+2(d 0 , d)> 3) ), and d® is determined 
from the state vector using the formulas by taking time derivatives of the motion equations 

in). 

After this substitution we get a 0 . 

A„ B n 


The correlation matrix P 0 = 
mined as follows: 

CTq, ( 7 a Up erg CJV 2u a j r 


B' c 0 


of the initial approximation errors a 0 is deter- 


^=2crg/r 


V= x 0 +y 0 +z 0 < 7 ^=a d +TV- max (u a , up) /A u~ d =2u d /r 

0(x y zj 

G 0 = 7-)—2— — (d Q ) is calculated using formulas following from (10) 
o(d, a, p) 


(26) 


1 

( a 2 - 0 0 \ 


/ a? 0 0\ 

d 

o 

© 

PCh 

0 e% 0 

R-o.i = 

0 a? 0 

a 

1 

0 0 P) 


0 0 u 2 . 

\ 0 / 

> 

o 

0 

o 

fS 

o 

o 

0 

o -■* 

a 

o 

II 

: G 0 R 01 G' 0 Bo — 0 


Here cr d , <j a , up - are the RMS values of typical errors of the measured parameters. 
Accounting for the system noise follows [6]. The matrix I\ is calculated in this way: 

110-^/(r 2 (Af efid ) 4 ) 0 0 

0 110-er 2 /(r 2 (Af ef , a) 4 ) 0 

0 0 110-^/(r 2 (At ef ^) 4 ; 


^ _ d(x, y, z) (A 

^ k ~ ~ 


d(d, a , (3 ) 


o o 
o r u i 


where 


- u d , a a , up - are the RMS values of typical values of the errors of the measured parameters, 

- Atgf d, Af e f ja , At e f t p - efficient memory of the algorithm for parameters d, a, f3 (param- 
eters of the algorithm selected on experimental basis). 



The cycle of processing the measurements d k , a k , /3 k (k=3,4,...) uses the formulas (23). 
They provide the following: using the previous parameters estimate a k ^ 1 and the calculated 
correlation matrix of its errors P fc _i, related to the time t k _ 1: and using the current scalar 

measurement u k , acquired for the time t k >t k _ i, calculate the updated estimate a fc and 
calculated correlation matrix of its errors P fc . First using (23) we update the previous 
estimate using u k =d k , then the estimate is additionally updated using u k =a k and finally 
the estimate is updated using u k =j3 k . The vectors h fc used in (23) are calculated using 
formulas resulting from partial differentiation of relationship (10). 


4. Non-recursive Method Based On the Guarantee Method 


The basic assumption of the Guarantee approach for non-linear problems is the small level 
of measurement errors and the known values of the upper limits of these errors. The term 
’’Guarantee” in this case means that the algorithm provides not only the calculated orbital 
parameters, but also the maximum possible errors of these parameters. The principals of 
the approach are as follows. 

Assume at the times t k (t k <t k+ p k — 1,2, ...n) we acquire measurements u k of certain 
functions h k { c) of the m- vector of parameters c (m<n), and the errors of the measurements 
5u k —u k —h k ( c) are limited from above by the constant 8 k . max . The estimate c n of parameters 
c obtained using the Guarantee approach and the vector 8c n;max of the maximum errors of 
the components of this estimate have the following geometrical interpretation. 

The limits of the measurement errors define in the m-dimensional space of parameters 
c = (Ci, c 2 , ..., c m ) a domain 

n 

11 n $i,max ^ h k (C) ^ U k 8 kmax ]■ (27) 

k= 1 


of possible values of c. We project this domain to the coordinate of the components of 
vector c and among the projected points for each axis find the most right c n . r and the most 
left c n i . They define the boundaries (maximum and minimum values) for the changes of 
each component of parameter c. In this case the estimate c n of parameter and maximum 
errors 8c n;rnax of this estimate are naturally defined as 


Cn 



5c 


n,max 



(28) 


If the measured parameters are linearly related to the determined parameter a, i.e. 
hi( a)=h'-a, where the components 6x1 of vector h* do not depend on a, this problem can 
be formulated as a standard problem of linear programming. Its solution is obtained in [22] , 
we will not describe the algorithm here. 
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If the measured parameters are linearly related to the determined parameter i.e. 
h k (c)=h' k -c, where components mxl of vector h fc do not depend on c, this prob- 
lem can be formulated as a standard promlcm of linear programming. The solution of the 
problem is known but the procedure is rather sophisticated and we will not present it here. 


For evaluation of the accuracy characteristics of the parameters obtaned using method it 
will be useful to consider a model case: we measure a scalar parameter c, with h — 1 and 
5k, max = 5max • In this one dimensional model 

n 

^maxj ^k~^~^max]} [ Hl&X Uk ^maxi Hlill ^k~^~^max ] (29) 

11 k k 

k= 1 

c n = \ ' (max u k + min u k ) 6cn :max =5 max -}--(maxu k - mmu k ) (30) 

k k k k 

where [a, b] - denotes an interval with left end a and right end b. 

The estimate c n from (30) has some important properties which to a certain extent are 
retained for the multi-dimensional case. 


1) The estimate (30) does not depend on 5 max , that means it is not sensitive to the accuracy 
of the knowledge of the characteristics of the errors - in this case to the accuracy of the 
value 5 max . 

2) If we know 5 max precisely the value 5c n _ max is a correct upper estimate for the error of the 
estimate c n , for any correlations of the errors of different measurements. Correctness in 
this case means the following. On one hand the true errors of the estimate a n for any 
n can not be greater than 5a n)rnax . On the other hand - the change of 5c n ^ max with 
growth of n correspond to the change in the correlation characteristics of the measurement 
errors. If the correlation interval of the errors is limited, the value 5c n rnax with the growth 
of n can be made indefinitely small. If this is not the case, i.e., a systematic error is 
present in the measurements, the value 5c n _ max for any n will be limited from below by 
the value of this error. 

3) For uncorrelated measurement errors the estimate c n with regard to accuracy is not inferior 
and sometimes can be essentially more accurate (!) than the estimate c n of the least squares 
technique, which in this case has the form 



k= 1 


(31) 


Thus, for example for the most frequent in practice uniform and triangular distributions of 
the measurement errors within the interval (—5 max ,5 max ) the estimates c n and c n are 
unbiased, and the r.rn.s. deviation of their errors and cr 5n for n>10 is calculated 
with a relative error not more than 10% using the following asymptotic formulas 
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for uniform distribution: 


1.4-5. 


n 


0.58-o maa; 



(32) 


- for triangular distribution: 


0,46-5. 




n 






0.41-5 


max 



(33) 


that follow from precise relationships, given in [5]. 

From (32) and (33) one can see that for a triangular distribution the accuracy of the esti- 
mates a n and a n is close and for the uniform one the estimate a n is about ~0 Ay/n times 
more accurate. Thus for the uniform distribution, not satisfying the regularity condition of 
Duget [5]"*", the estimate a n turns to be super-efficient. 

These properties make the Guarantee approach more advantageous compared to the least 
squares method and make it rather attractive for use in applications characteristic for the 
operation of different radars. This is especially characteristic for rather accurate and stable 
in performance radars where the abnormal measurements (if any) can be easily selected dur- 
ing preliminary processing of the single measurements and do not enter the process of orbit 
determination. 


Let us treat the considered task in more detail. 

Assume that the radar at the times t k (k = 1,2 measures the range d, azimuth 
oi and elevation angle (3 in the local spherical coordinate frame and the errors of the 
measurements do not exceed respectively the values 5d,max, da, max, 8p,max- The task is to 
determine at the time t=0.5-(fi+f n ) the six-dimensional vector of orbital parameters c = 
( d , a, /3, d, a, (3) and its maximum errors in this coordinate frame. 

The following algorithm for solving the task is suggested. 

Divide all the measurements into nj 2 groups*. The first group includes the measurements 
acquired at L and t 05n , the second group - the measurements with t 2 and t 0 , 5n+1 , etc. 
For each k - th group (k = 1,2, ,.,0.5-n) from two position vectors u k = (d k , a k , /3 k ) and 
u 0 . 5 n+fc = (d 0 . 5 n+fc, a 0 . 5 n+fc, P 0.5 n+k) in the local spherical coordinate frame we determine the 
six dimensional vector of orbital parameters c fc = (d, a, (3, d, a, $) k in the same coordinate 
frame for the time i. 

Instead of the labor consuming operation of determination of the domain (27) and its 
projections on the coordinate axes of the phase vector we suggest to project on these axes the 

Tn this case the break points of the distribution depend on the estimated parameter. 

*We assume that n is even. If n is odd we consider that for the time t n we have two identical 
measurements. 
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domain D fc , corresponding to the k-th group and then for each axis search for the intersections 
of this projections. 

The domain D fc is approximated by a six-dimensional parallelepiped with the center in 
the point c k and 64 apexes determined by the formula 


c, ± S, 


d,max \ 


ji ± 5 , 


a, max 


■h 


± Sr-, 


(3, max' J 3 ^ djnax ‘ Ji 0 arnax - J 5 d(3,max' J6 


± S n 


± Sr-, 


(34) 


where j 1; j 2 , j 3 , j 4 , j 5 , j 6 - are the lines of the matrix of partial derivatives J = 
(9(C^ ) 

- of the functional transformation c fc =f fc (u fe , u 05n+fc ) in the point c^. The 

o( u fe , u 0 . 5n+fc ) 

boundary projections c k)l and c k , r of these 64 points of the six dimensional phase space to 
each of its axes determines the boundaries of the vector interval (c M , c fcjr ) of the possible 
values of all six orbital parameters determined by the k-th group of measurements. 

n 

After determination of the common part (c ; , c, r ) = f) (c fc)i , c kjr ) of these intervals for all 

k = 1 

groups of measurements using eq. (28) we determine the estimate of the orbital parameters 
and its maximum errors. 

The major element of the suggested algorithm is the procedure c fc =f fe (u fc , u 0 . 5n+fc ) for 
orbit determination on the basis of two positions. This procedure is used for determination 
of the state vector c fc , and the matrix of derivatives J. A rather simple and efficient 
algorithm for solving this task for Kepler’s approximation is given in [7] (algorithms 31 and 
51). The algorithm uses the measured position vectors r = (A", Y, Z) of the satellite 
in a fixed coordinate frame for Kepler’s orbit. They are obtained from the measurements 
u fc = ( d k , a k , (3 k ) in the following way: 

- For the fixed coordinate frame X, Y, Z we select the system whose axes for the time 
t = 0.5(tj+t n ) coincide with the axes of the Greenwich system. The measurements 
u fc = ( d k , a k , (3 k ) are transferred to r k = (X k , Y k , Z k ). 

- From all the measurements r k we substraet such values of Ar k , for which r k =r k — Ar k 
follow (up to the errors of the measurements) Keplerian orbit, whose parameters for the 
time t = 0.5-(ti+t n ) coincide with the parameters of real (perturbed) orbit. To obtain 
Ar k we solve the equations of motion using Euler’s method 


Ar = -fi - at + 3-/i-(Ar, r)-r — 3-jk-R 2 - J 2 -(0.5-r + zuj 0 — 2.5-z 2 -r ) (35) 

for the deviations of Ar, Ar from the Keplerian orbit for Ar ( t ) = Ar ( t ) = 0 and r = r k . 
Here r = \f (r, r) ? = r/r Ar = Ar /r z = z/r fx = fi/r 2 R=R e /r. 
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5. Conclusions 


A new method for determining the orbit of a space object from a single pass through a 
radar has been presented. This method, called the Guarantee Method, is based on linear 
programming and offers advantages over the standard least squares and Extended Kalman 
filter. Under certain conditions this method provides a much more accurate state estimate 
than the least squares and Kalman filter. It also provides limits on the state estimate errors 
when the measurement errors have a uniform or triangular distribution. The amount of 
computation is significantly more for the Guarantee Method than the other methods, but this 
is not a problem with modern computers. A comparative analysis with the other methods has 
been presented. 
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